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Abstract:  Simulated  bidirectional  reflectance  distribution  functions 
(BRDF)  were  compared  with  measurements  made  just  beneath  the  water’s 
surface.  In  Case  1  water,  the  set  of  simulations  that  varied  the  particle 
scattering  phase  function  depending  on  chlorophyll  concentration  agreed 
more  closely  wtth  the  data  than  other  models.  In  Case  II  water,  however,  the 
simulations  using  fixed  phase  functions  agreed  well  with  the  data  and  were 
nearly  indistinguishable  from  each  other,  on  average.  The  results  suggest 
that  BRDF  corrections  in  Case  II  water  arc  feasible  using  single,  average, 
particle  scattering  phase  functions,  but  that  the  existing  approach  using 
variable  particle  scattering  phase  functions  is  still  warranted  in  Case  1  water. 
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OCIS  codes:  (010.4450)  Oceanic  optics,  (010.4588)  Oceanic  scattering;  (010.5620)  Radiative 
transfer,  (280.4991)  Passive  remote  sensing. 
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1.  Introduction 

The  spectral  radiance  exiling  a  natural  water  body  (the  water-leaving  radiance,  L»)  is 
dependent  on  the  solar  zenith  angle  and  the  viewing  direction,  as  well  as  the  optical  properties 
of  the  water  and  its  constituents.  When  this  radiance  is  imaged  with  sensors  on  earth -orbiting 
satellites  (ocean  color  remote  sensing),  every  pixel  in  the  image  is  associated  with  a  different 
illumination-viewing  geometry.  Most  algorithms  that  are  used  to  associate  water-leaving 
radiance  with  geophysical  parameters,  such  as  the  concentration  of  chlorophyll  a ,  have  been 
developed  based  on  in  situ  measurements  of  upwelling  spcetral  radiance  propagating  toward 
the  zenith  [1],  i.e.,  at  a  single  viewing  angle.  Thus,  it  is  neeessary  to  be  able  to  relate  the 
water-leaving  radiance  at  any  sun-viewing  direction  to  the  nadir  view.  In  addition,  today  there 
are  many  earth-orbiting  sensors  in  operation,  and  each  views  a  given  location  under  different 
illumination-viewing  geometries.  Therefore,  in  order  to  compare  or  merge  the  imagery  from 
different  sensors,  it  is  neeessary  to  relate  each  to  a  common  geometry  [2]. 

Gordon  and  Clark  [3]  first  defined  the  normalized  water-leaving  radiance  (nL*)  as  the 
value  would  have  if  the  sun  were  at  the  zenith,  the  atmosphere  absent,  and  the  sensor 
looking  at  nadir.  This  estimate  was  obtained  by  assuming  that  the  angular  distribution  of 
upwelling  radiance  just  beneath  the  water’s  surface  (the  bidirectional  reflectance  distribution 
function,  BRDF)  was  independent  of  the  viewing  direction,  i.e.,  they  assumed  a  completely 
diffuse  BRDF.  Previous  efforts  to  model  /  understand  the  actual  BRDF  [4-10)  have  produced 


#160488  -  $15.00  USD  Received  23  Dec  201 1;  revised  5  Mar  2012;  accepted  5  Mar  2012;  published  20  Mar  2012 
(C)  2012  OSA  26  March  2012 /Vol.  20,  No.  7  /  OPTICS  EXPRESS  7631 


the  general  insight  that  the  variation  can  be  modeled,  at  least  in  Case  I  waters  (those  for  which 
the  optical  properties  eovary  with  the  chlorophyll  concentration  [  1  ], ). 

The  Morel  et  al.  [9)  model  (referred  to  below  as  MAG2002)  is  currently  the  standard 
model  used  to  account  for  the  angular  variation  of  the  upwelling  radiance  in  satellite 
processing.  It  works  well  in  Case  1  water  [11],  but  is  not  suitable  for  Case  II  water  for  at  least 
three  reasons,  hirst,  the  tables  end  at  a  maximum  chlorophyll  concentration  of  10  pg  \  \  but 
chlorophyll  can  exceed  this  in  places.  Second,  the  MAG2002  BRDF  tables  were  developed 
from  a  radiative  transfer  (RT)  model  that  used  scattering  particle  phase  functions  that  assumed 
that  all  the  particles  arc  biogenic,  and  thereby  are  applicable  only  to  Case  1  waters.  Third,  it 
requires  that  the  absorption  contribution  from  other  optically  active  agents  (e.g.  yellow 
substance)  co-vary  with  concentration  of  chlorophyll  a%  which  is  often  invalid  for  Case  II 
waters. 

Radiative  transfer  computations  require  the  absorption  coefficient  (a)  and  the  volume 
scattering  function  where  6  is  the  scattering  angle),  or  equivalently,  the  absorption 

coefficient,  the  scattering  coefficient  (/>),  and  the  scattering  phase  function,  =  0(0)/ b). 
Of  these,  measurements  of  p(0)  or  the  phase  function  have  been  reported  in  only  a  few 
studies,  most  of  which  were  reviewed  by  [12]  and  references  therein.  Recently,  using  novel 
instrumentation,  Sullivan  and  Twardowski  [13]  found  remarkable  consistency  in  the  shape  of 
the  backward  portion  of  /?  (i.c.  for  90°  <  =  (?<=  180°)  measured  in  situ  around  the  world. 
This  is  a  promising  development,  because  a  “universal”  particulate  scattering  phase  function 
would  simplify  radiative  transfer  modeling  in  coastal  waters.  In  particular,  combining  this 
phase  function  with  a  and  bb  (the  baekseattering  coefficient),  which  are  accessible  from  the 
remote  sensing  signal  [14],  would  enable  the  BRDF  to  be  determined. 

One  objective  of  this  study  was  to  assess  the  BRDF  corrections  from  a  recently  published 
algorithm  [10],  which  may  apply  in  cither  Case  I  or  Case  II  water.  Lee  et  al.  [10]  (referred  to 
below  as  Lec201 1)  used  in  situ  measurements  from  just  3  locations  to  validate  their  model; 
here  we  used  a  much  larger  data  set  across  a  wide  variety  of  inherent  optical  properties. 

A  second  objective  of  this  study  was  to  assess  the  effects  of  using  different  phase 
functions  when  modeling  the  BRDF.  MAG2002  and  Lce201 1  used  different  phase  functions. 
MAG2002  used  a  phase  function  that  varied  as  a  function  of  chlorophyll  concentration.  In 
contrast,  the  phase  function  used  by  Lee20l  1  did  not  vary  systematically  as  a  function  of  the 
BRDF  look-up-tablc  input  parameters.  However,  the  synthetic  data  set  used  to  generate  the 
Lee201 1  BRDF  look-up-table  had  been  computed  with  a  randomly  varying  linear  mixture  of 
two  phase  functions  (one  for  mineral  particles,  one  for  phytoplankton).  One  question  posed  in 
Lee201 1  but  not  addressed  with  experimental  data  was  how  critical  the  choice  of  phase 
function  is  for  the  purposes  of  a  BRDF  correction.  We  investigated  that  question  here  by 
comparison  of  the  Lec201 1  BRDF  correction  with  those  provided  by  MAG2002  as  well  as 
our  own  corrections  generated  with  a  radiative  transfer  model  and  the  Sullivan  and 
Twardowski  [13]  and  Petzold  [15]  turbid-water  phase  functions. 

The  overall  goal  of  this  work  was  to  assess  potential  BRDF  corrections  applicable  in  Case 
II  water.  As  mentioned  above,  two  specific  objectives  were  (a)  to  validate  Lec201 1  using  a 
large  data  set  of  in  situ  hemispherical  radiance  measurements  and  (b)  to  test  the  sensitivity  of 
results  to  the  shape  of  the  phase  function.  For  comparison,  wc  also  included  data  from  Case  I 
sites  and  corrections  from  MAG2002  in  the  evaluation.  Although  the  focus  was  on  Case  II 
water  and  the  new  model  Lec201 1,  the  older  model  MAG2002  provided  a  useful  benchmark, 
since  it  had  been  previously  validated  in  Case  1  conditions  using  methods  similar  to  those 
presented  here  [II]. 

2.  Methods 

The  approach  was  to  compare  measured  hemispherical  upwelling  radiance  distributions  with 
model  simulations.  Each  comparison  consisted  of  a  set  of  average  images  measured  at  412, 
436,  486,  526,  and  548  nm  with  the  NuRADS  system  [16],  and  sets  of  results  from  each  of 
five  numerical  models.  The  first  set  of  model  output  was  interpolated  from  the  tables 
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described  by  MAG2002.  The  second  set  of  model  output  was  interpolated  from  the  tables 
described  hy  Lee 20 1 1 .  Three  additional  sets  of  model  output  were  generated  using  our  own 
RT  code  (H.R.  Gordon,  unpublished)  with  different  fi(Q)  parameterization,  Data  were 
availahle  from  six  field  campaigns:  Chesapeake  Bay,  USA  in  2004,  the  BIOSOPE  cruise  in 
2004  [17],  Monterey  Bay,  CA,  USA,  in  October  2006,  the  Ligurian  Sea  in  both  October  2008 
(LSCV’08)  and  March  2009  (BP’09),  and  the  New  York  Bight  in  2009. 

2.  /.  NuRADS  image  processing 

The  NuRADS  system  [16]  is  a  compact  (30  cm  diameter,  30  cm  length),  multispectral  camera 
that  images  the  upwelling  light  field  in  six  narrow  (-10  nm  FWHM)  spectral  bands  centered 
at  412,  436,  486,  526,  548  and  615  nm.  This  study  did  not  use  data  from  the  615  nm  channel 
due  to  the  large  influence  of  instrument  self-shading  at  this  wavelength. 

NuRADS  was  deployed  on  the  surface  with  a  nadir  view  direction  thereby  collecting 
upwelling  radiance  -30  cm  below  the  surface.  The  camera  sequentially  acquired  images  in 
each  hand  with  the  use  of  a  rotating  filter  changer.  Typical  exposure  times  were  less  than  one 
second.  Acquiring  a  set  of  images  from  all  six  wavelengths  took  about  two  minutes,  however, 
due  to  the  delay  reading  the  data  from  the  CCD  in  between  exposures  for  each  spectral  hand. 
A  typical  deployment  lasted  from  one  to  several  hours,  enabling  multiple  acquisitions  over  a 
range  of  solar  zenith  angles. 

Reduction  of  the  raw  NuRADS  images  consisted  of  applying  calibration  factors  (see  [16]) 
and  averaging  images  in  both  space  and  time  to  reduce  environmental  noise.  Before 
averaging,  every  image  was  inspected  to  find  the  anti-solar  point,  to  correct  the  geometry  of 
the  image,  and  to  check  for  obstructions  in  the  field  of  view  such  as  fish,  the  power/data 
cable,  the  side  of  the  ship,  or  other  instruments.  Images  were  then  averaged  in  10-minute  bins, 
excluding  those  that  had  been  flagged  as  unacceptable  in  the  inspection  stage.  The  symmetry 
of  the  images  about  the  principal  plane  was  exploited  to  further  average  both  halves  of  each 
image.  In  addition,  spatial  binning  over  3x3  pixel  windows  was  performed  to  produce  final 
average  images  at  1  x  I  degree  resolution.  Therefore,  each  pixel  in  the  processed  image  over  a 
10-minute  window  for  a  given  hand  could  be  an  average  of  up  to  90  raw  pixels  (5  images  x  2 
image  halves  x  9  pixel  window).  The  mean  and  coefficient  of  variation  (CV  *  standard 
deviation  divided  by  the  mean)  were  computed  for  each  pixel  in  the  processed  image  from  the 
up  to  90  raw  pixels  in  the  original  images. 

2.2.  I  OP  data  processing 

Inherent  optical  properties  (lOPs)  measured  during  the  field  experiments  along  with  solar 
zenith  angle,  calculated  for  the  time  of  each  NuRADS  image,  were  used  to  index  the  look  up 
tables  (Sections  2.3,  2.5)  and  parameterize  the  RT  model  (Section  2.6).  The  solar  zenith  angle 
((?,),  defined  as  the  angle  between  the  vector  to  the  sun  and  local  vertical,  was  computed  based 
on  the  time,  latitude,  and  longitude  of  each  NuRADS  image.  Absorption  ( aw )  and  scattering 
(bH)  coefficients  of  seawater  were  interpolated  to  NuRADS  spectral  bands  from  Table  l.l  in 
[18].  The  absorption  coefficient  of  dissolved  and  particulate  constituents  (an)  and  the  particle 
total  scattering  ( bp )  and  backscattering  ( b ^)  coefficients  were  measured  in  situ ,  depth 
weighted,  and  interpolated  as  necessary  to  the  NuRADS  spectral  bands  as  follows. 

Vertical  profiles  of  the  particulate  and  dissolved  absorption,  apg(A,  z),  and  particle 
scattering  bp(k,  z)  coefficients  were  measured  at  9  wavelengths  (x)  using  WET  Labs  ac-9 
instruments  during  the  Chesapeake,  BIOSOPE,  LSCV’08,  BP’09,  and  New  York  Bight 
experiments  and  at  83  wavelengths  using  a  WET  Labs  ac- s  in  Monterey  Bay.  Particle 
backscattering  coefficients  6^(A,  z)  were  measured  with  ECO-BB3  instruments  during 
BIOSOPE,  Monterey  Bay,  LSCV’08,  and  New  York  Bight  cruises  and  with  a  HOBI  Labs 
Hydroscat-6  instrument  for  BP’09.  No  b ^  measurements  were  made  during  the  Chesapeake 
cruise.  Quality  control  and  calibrations  were  applied  to  these  data  sets  by  the  investigators 
who  acquired  them  [  1 9,  20]. 

Processing  the  calibrated  apg(X,  z),  bp(X,  z),  and  bhp(X ,  z)  followed  the  same  four  steps  for 
each  cast  in  each  of  the  data  sets.  First,  we  computed  depth-weighted  values  b^ rtghteJ^mn) 
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and  «/)  using  Eq.  (I)  for  the  shortest  wavelength  data  (A*,,*)  from  the 

baekseatteri ng  sensor  and  the  closest  corresponding  ac-9  or  ac-s  band  (kctoSf)  in  each 
experiment.  The  k ^  /  A</Mr  bands  were  470  /  488  nm,  450  /  45 1  nm,  462  /  488  nm,  442  /  440 
nm,  and  462  /  488  nm  in  the  BIOSOPE,  Monterey  Bay,  LSCV'08,  BP*09,  and  New  York 
Bight  data  sets,  respectively. 


(I) 


In  Eq.  (1),  v(A,  z)  is  the  depth-dependent  variable  being  averaged  (either  apg(k%  z),  bp(k ,  z), 
or  bty(k,  z)),  K*ixhir*i(k)  is  the  depth- weighted  value  of  that  variable,  a, (A,  z)  and  b^k,  z)  are  the 
total  depth-dependent  absorption  (akk,  z)  =  tf^A,  z)  +  tf*(A))  and  baekseatteri  ng  (6*(A,  z)  = 
btyi A,  z)  +  fym(A))  eoeffieients,  respectively,  and  Zmax  is  the  maximum  depth  of  the  east. 

Second,  we  used  the  baekseatteri  ng  ratio,  Bp  *  /  and  the 

assumption  that  Bp  did  not  vary  with  wavelength  [19,  21]  to  compute  6^(A,  z)  =  BJbp(k,  z)  for 
all  of  the  ac-9  /  ac-s  spectral  bands.  Third,  we  computed  depth-weighted  absorption  apg 
*w*ai*AA),  scattering  V»w,w(A),  and  baekseattering  b^^M)  eoeffieients  using  Eq.  (1)  for 
all  of  the  ac- 9  /  ac-s  spectral  bands.  Finally,  the  depth-weighted  values  were  interpolated  to 
the  NuRADS  wavelengths.  Linear  interpolation  was  used  for  an^Tij(^.  A  power  law  of  the 
form  ak  was  used  to  interpolate  b^^  and  Note,  all  further  references  to  an,  b ^ 

or  bfy  imply  the  depth-weighted  value  of  those  variables  interpolated  to  a  particular  NuRADS 
spectral  band. 

2.3.  Chlorophyll  processing 

Total  chlorophyll  concentration  [Chi]  (pg  1  !)  represents  the  sum  of  the  concentrations  of  the 
following  suite  of  pigments:  chlorophyll  a,  divinyl  chlorophyll  a,  ehlorophyllid  a,  and 
chlorophyll  a  allomers  and  epimers.  Calculation  of  [Chi]  was  determined  via  High- 
Performance  Liquid  Chromatography  (HPLC). 

Calculation  of  [Chi]  from  the  BIOSOPE  cruise  used  a  slightly  modified  version  (see  [22]) 
of  the  method  initiated  by  [23].  Data  for  the  other  cruises  were  processed  following  NASA 
protocols  [24].  All  of  the  Chesapeake  water  samples  were  acquired  at  0.5  m  depth.  Water 
samples  for  the  other  cruises  were  acquired  within  2  m  of  the  surface  and,  for  some  stations, 
at  depths  as  great  as  20  m.  Stations  for  which  a  depth  east  was  available  were  optically 
weighted  using  Eq.  (1)  to  generate  a  single  [Chi]  value. 

2.4.  MAG2002  table  interpolation 

MAG2002  provided  look-up  tables  for //  Q  as  a  function  of  wavelength,  solar  zenith  angle, 
chlorophyll  concentration,  and  view  angle.  One  set  of  MAG2002  look-up  tables  included  the 
effects  of  Raman  scattering  whereas  a  second  set  of  tables  did  not.  The  results  here  used  the 
tables  that  did  include  Raman  scattering.  The  dimensionless  parameter/ is  proportional  to  the 
irradianee  reflectance;  /  =  (a,  /  bN)  *  (E„  /  Ed ),  where  E„  and  Ed  are  the  upwelling  and 
downwelling  irradianees  just  below  the  water  surface,  respectively.  The  bidirectional  function 
Q  is  defined  as  the  ratio  of  the  upwelling  irradianee  to  radiance;  Q  =  (£*  /  L J,  where  Ly  is  the 
upwelling  radiance  just  below  the  water  surface. 

Total  chlorophyll  concentration  and  solar  zenith  angle,  calculated  for  the  time  of  each 
NuRADS  image,  were  used  to  retrieve  f  /  Q  at  regularly  spaced  in-water  view  ( 0V )  and 
azimuth  (0)  angles  via  interpolation  of  the  MAG2002  look-up  tables.  Azimuth  is  defined 
relative  to  the  principal  plane,  the  plane  containing  the  anti-solar  point,  local  vertical,  and  the 
sun.  The  tables  were  not  interpolated  spectrally  because  they  had  been  computed  at 
wavelengths  close  to  the  NuRADS  bands  (412.5,  442.5,  490,  510,  560  nm).  Normalizing/ /  Q 
at  each  view  angle  by  the  value  at  nadir  ( 0V  =  0,  0=0),  as  described  in  Section  2.7,  produced  a 
normalized  radiance  distribution,  LJ,0^  0)  /  LJ, 0,  0),  for  comparison  with  the  normalized 
NuRADS  images. 


#160488  -  SI 5  00  USD  Received  23  Dec  201 1;  revised  5  Mar  2012;  accepted  5  Mar  2012;  published  20  Mar  2012 
(C)  2012  OSA  26  March  2012  /Vol.  20,  No.  7  /  OPTICS  EXPRESS  7634 


2.5.  LeclOll  table  interpolation 

Lcc20l  I  provided  look-up  tables  for  dimensionless  parameters  G0  and  G/,  which  relate 
remote  sensing  reflectance  to  the  water  and  particle  plus  dissolved  matter  absorption  and 
backscattering  coefficients  via: 


( 0. .  0m  ) -  (  g;  (0, .  .* ) + g;  (O, .  o„  »  ^ ^ 


(2) 


Jl 

K 


where  k  =  a,  +  and  G0  and  G;  arc  empirical  coefficients  determined  separately  for 
water,  G*  and  particles,  G'1,  with  dependence  on  solar  zenith  angle,  above-water  view  angle 
(0ltJ),  and  azimuth.  Tables  of  coefficients  G  were  defined  for  0S  <  =  75°.  Remote  sensing 
reflectance,  Rn ,  is  defined  as  the  ratio  of  the  water  leaving  radiance,  to  the  downwclling 
irradiance,  measured  just  above  the  surface. 

Normalizing  Rn  at  each  view  angle  by  the  value  at  nadir  ( 0 w  =  0,  <f>  =  0)  produced  a 
normalized  above-water  radiance  distribution,  L*(0i%  0 w,  $  /  0,  0).  Note  that  is 

canceled  when  Rn  is  normalized  in  this  way.  Correction  for  the  reflection-transmission 
properties  of  the  air-sea  interface  was  necessary  to  compare  normalized  with  normalized 
which  is  the  quantity  measured  by  NuR  ADS: 


=  Lw(0„d„,+)  9UQ) 

L.(0,. 0.0)  LJ0„ 0,0)  tt<0 

The  dimensionless  factor  91  in  principle  depends  on  both  Q,  and  0TO,  but,  as  shown  in  [25], 
9?  can  be  approximated  by: 


*(0 


0.957  »  7)  (0,J 
0.985  *m2 


(4) 


to  within  a  few  percent  for  0X*  <  60 }  and  to  better  than  10%  for  0iXt  =  70°.  In  Eq.  (4),  T/,0^)  is 
the  Fresnel  transmittance  for  a  flat  air-water  interface  and  in  is  the  refractive  index  for  water, 
taken  to  be  1.341 . 

Solar  zenith  angle,  a apv  bh„  and  b^  measured  and  processed  as  described  in  Section  2.2 
were  used  to  retrieve  R„  at  regularly  spaced  angles  ( 0 ,a,  <f>)  via  interpolation  of  the  Lce201 1 
look-up  tables  (i.e.  Rn  was  derived  from  lOPs  and  the  tables,  not  measured  directly  in  the 
field).  Rn(0% „,  <f>)  was  converted  to  normalized  Lu(0„  <f>)  for  comparison  with  the  NuRADS 
images.  As  mentioned  in  Section  2.2,  no  b ^  measurements  were  made  dunng  the  Chesapeake 
cruise,  yet  b ^  was  required  for  F^q.  (2).  To  address  this  limitation,  we  assumed  Bp  =  2%,  a 
common  value  for  coastal  waters  [21,  26],  and  computed  b^  =  Bp*bp.  To  check  the  validity  of 
this  assumption,  we  ran  the  same  analysis  using  Bp  =  1%  and  3%  and  computed  normalized  l * 
from  Eqs.  (2-4)  using  the  measured  bp  and  Bp=  \  %  and  3%.  The  normalized  L*  values  for  Bp 
=  1%  and  3%  agreed  within  3.8%  for  95%  of  the  Chesapeake  data  set  and  within  5.6%  for 
99%  of  the  Chesapeake  data  set. 

2.6.  Our  RT  model  parameterization 

Our  RT  code  required  the  following  input  parameters:  the  lOPs  and  solar  zenith  angle 
described  in  Section  2.2,  the  Rayleigh  optical  depth  of  the  atmosphere,  computed  from  [27], 
and  a  particle  volume  scattering  function,  flp(0).  Varying  the  choice  and  /  or  normalization  of 
Pfj[0)  generated  three  sets  of  model  output  for  comparison  with  each  NuRADS  image. 

Method  I  used  the  Petzold  [15]  turbid-water  fip(0)  as  normalized  in  [28]  to  produce  the 

Petzold  turbid-water  phase  function,  P  \0) .  Note  that  this  pp(0)  satisfies  the  normalization: 
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(5) 


2n\Pp(O)s\n{9)d0=\ 

0 

To  recreate  the  proper  p^O)  for  input  to  the  RT  code  under  method  1  we  multiplied  p  (0) 

by  the  depth-weighted,  spectra  I  ly-interpolatcd  bp.  Note  that,  due  to  the  acceptance  angle  of  the 
ae-9  (0.93  degrees),  our  depth-weighted  bp  underestimates  the  value  that  would  be  measured 
by  a  perfect  instrument  for  the  Pet/old  phase  function  by  20-30#  [15]. 

Method  2  used  the  Sullivan  and  Twardowski  [13]  average  phase  function.  The  MASCOT 
instrument  used  in  [13]  collected  data  over  the  range  from  10  to  170  degrees,  but  Sullivan  and 
Twardowski  extrapolated  the  average  phase  function  in  the  backwards  direction  from  170  to 
180  degrees  [13].  Extrapolation  of  the  phase  function  in  the  forward  direction  is  more 
difficult;  therefore  we  performed  some  initial  experiments  to  test  the  importance  of  the  exact 
shape  of  the  forward  direction  on  the  BRDF.  These  tests  consisted  of  a  set  of  model  runs  at 
412  nm,  over  a  range  of  apv  and  bp  from  0.01  to  1 .0  m  1  and  solar  zenith  angles  0,  30,  and  60 
degrees,  in  which  we  truncated  the  Petzold  phase  function  at  2,  5,  and  10  degrees.  The 
normalized  upwclling  radiance  distributions  created  from  the  phase  function  truncated  at 
either  5  or  10  degrees  were  always  within  +/-  3#  of  the  corresponding  result  using  the  phase 
function  truncated  at  2  degrees.  Therefore,  we  did  not  extrapolate  the  Sullivan  and 
Twardowski  phase  function  in  the  forward  direction,  instead  using  it  truncated  at  10  degrees 
(Fig.  I). 

The  observation  that  truncating  the  phase  function  does  not  affect  the  calculated  upwclling 
radiance  significantly  is  consistent  with  Mobley  et  al.  [29],  who  concluded  that  the  exact 
shape  of  the  phase  function  in  the  forward  direction  was  not  critical  as  long  as  bhP  was  correct. 
Furthermore,  Gordon  [30]  showed  that  P(0)  truncated  at  scattering  angles  <  15  degrees  (such 
as  the  one  used  here)  did  not  alter  calculations  of  irradianees  by  more  than  a  few  percent. 

Note  that  Sullivan  and  Twardowski  [13]  used  a  normalization  in  the  backward  direction 
only,  i.e.  their  measured  pf/0)  were  normalized  by  b^  and  not  by  bp  as  had  been  done  in  [28]. 

Therefore,  the  pp(6)  for  method  2  satisfies  the  normalization: 

M 

In  J  Pf(d)%\nW)dd=\  (6) 

Mil 

To  recreate  the  proper  fi^O)  for  input  to  the  RT  code  under  method  2  we  multiplied  p  (0) 

by  the  depth-weighted,  speetrally-interpolated  bhp. 

Method  3  also  used  the  Petzold  [15]  turbid- water /?(/?),  but  rather  than  scaling  it  by  br  as 

in  method  I,  it  was  normalized  in  the  backward  direction,  as  in  method  2.  Thus,  /?,(#)  for 

method  3  satisfied  Eq.  (6)  and  was  multiplied  by  the  depth -weighted,  speetrally-interpolated 
bty  to  produce  the  P(0)  input  to  the  RT  calculations.  The  p(0)  used  for  method  3  has  the  same 
shape  as  p(0 )  used  for  method  I ,  but,  like  method  2,  forces  b^  in  the  model  to  agree  with  the 
measurements. 
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Fig.  I.  Compariwo  of  Sullivan-Twardowski  and  Petzold  turbid  water  phase  functions  used  in 
methods  1-3.  The  phase  function  for  method  1  U  normalized  over  all  angles  (Eq.  (5))  and 
multiplied  by  b ,  to  generate  PJO)  in  the  RT  model.  The  phase  functions  for  methods  2  and  3 
arc  normalized  in  the  backward  direction  (Eq.  (6))  and  multiplied  by  b*  to  generate  fi/O)  in  the 
RT  model.  Inset  xhow%  details  of  methods  2  and  3  over  the  backscattering  directions. 

2.7.  Evaluation  of  the  model-data  difference  computation 

Each  of  the  radiance  distributions  was  normalized  by  its  value  at  nadir.  Comparison  between 
the  modeled  and  measured  normalized  radiance  distributions  was  performed  by  computing  the 
model-data  difference  at  every  5  degrees  in  nadir  from  5  to  45  degrees  and  every  15  degrees 
in  azimuth  from  0  to  180  degrees.  Thus,  the  difference,  D,  between  each  model  output, 
and  the  corresponding  average  NuRADS  image,  L */>,  was  computed  at  1 17  angles  based  on 
Eq.  (7). 


D(0vj) 


A.* (*,.*> 

A,*  (0.0)  4/>(0.0) 


(7) 


The  NuRADS  shadow  affected  some  of  the  view  angles,  which  were  subsequently 
discarded  based  on  an  estimate  of  the  shadow  contamination.  Instrument  self-shading  was 
estimated  using  the  Gordon  and  Ding  [3IJ  model  modified  to  incorporate  the  three- 
dimensional  shape  of  NuRADS  rather  than  the  flat  circular  disk  used  in  the  original 
calculations.  For  each  view  angle,  the  path  length  through  the  instrument  shadow  (ds)  was 
estimated  and  used  with  the  total  absorption  (at)  to  compute  the  transmission  ( T)  through  the 
shadow:  T  =  cxp(-arfs).  For  most  of  the  data,  angles  for  which  the  transmission  was  less  than 
95%  were  discarded  from  the  comparison.  The  total  absorption  coefficient  in  the  Chesapeake 
Bay  data  set,  however,  was  so  high  (mean  +/-  std.  dev.  =  3.1  +/-  1.3  m  *)  that  almost  the 
entire  data  set  was  flagged  as  shadow.  Therefore,  the  shadow  threshold  for  the  Chesapeake 
data  set  was  reduced  to  90%  transmission. 

After  omitting  points  thought  to  be  contaminated  by  instrument  shadow,  the  distributions 
of  the  remaining  D(0„  $)  were  plotted  as  grouped  by  chlorophyll,  and  0V  These  plots  were 
constructed  using  all  images  in  the  data  set  and  also  after  dividing  the  images  into  “Case  P 
and  “Case  11”  subsets  based  on  the  visual  relationship  between  aPK  or  bp  and  chlorophyll,  as 
described  below. 


3.  Results 

3.  /.  Case  I  and  Case  II  divisions  within  the  data 

The  optical  properties  of  Case  1  waters  are  determined  by  the  chlorophyll  concentration  [  1 1. 
Our  data  set  was  divided  into  Case  1  and  Case  11  subsets,  therefore,  by  visual  inspection  of 
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plots  of  inherent  optical  properties  as  a  function  of  chlorophyll  concentration  (Fig.  2).  The 
IOP  model  of  [6]  (their  Eq.  (7-10)  was  used  as  a  reference  during  this  process.  Although  the 
division  between  Case  1  and  11  was  done  visually,  in  effeet  the  criteria  for  classification  as  a 
Case  II  site  were  (a)  [Chi]  >  1 0  pg  1~\  or  (b)  bp  >  I  m-1,  or  (c)  an  >  0.2  m  *. 

The  entire  data  set  contained  1474  averaged  images,  80%  of  which  were  acquired  in  Case 
1  conditions.  The  example  shown  (Fig.  2)  is  for  526  nm,  which  had  319  averaged  NuRADS 
images.  Typically,  many  images  were  acquired  at  each  IOP  sampling  station,  which  is  why 
only  1 12  combinations  of  apx  and  bp  were  plotted  in  Fig.  2. 

3.2.  Mode l -data  comparison 


The  number  of  average  NuRADS  images  processed  and  matched  with  the  corresponding  runs 
from  our  RT  model  was  265,  284,  290,  319,  and  316  in  each  of  the  five  NuRADS  bands  (412, 
436,  486,  526,  and  548  nm,  respectively).  Of  these,  only  233,  254,  257,  284,  and  279  in  the 
respective  bands  were  matched  with  the  MAG2002  tables  because  the  remainder,  all  of  which 
were  from  cither  the  Chesapeake  or  Monterey  Bay  data  sets,  had  total  chlorophyll  values 
greater  than  10  pg  / 1  (the  upper  limit  available  from  MAG2002).  Also,  only  263,  279,  284, 
306,  and  305  images  in  the  respective  bands  were  matched  with  the  Lcc20l  1  tables  because 
the  remainder  were  acquired  with  solar  zenith  angles  >  75°. 
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Fig.  2.  Depth-  weigh  led  bf  (lop)  and  an  (bottom)  at  526  nm  plotted  against  local  chlorophyll 
concentration.  Left  plots  on  both  the  lop  and  bottom  .show  details  near  the  origin.  Solid  lines 
plot  lhc  IOP  modet  used  as  a  reference  for  Case  l  16].  Solid  dots  and  open  circles  represent 
conditions  designated  as  Case  l  and  Case  11,  respectively 


For  Case  I  conditions,  the  correlation  between  the  MAG2002  modeled  LJ,0„  <f>)  t  Ly( 0,  0) 
and  the  NuRADS  data  was  similar  to  that  between  the  Lec20l  1  model  and  the  NuRADS  data 
(Fig.  3A,  B).  The  Lce20l  I  model  more  closely  agreed  with  the  NuRADS  data  than 
MAG2002  in  Case  II,  however  (Fig.  3C,  D).  For  these  Case  I  data  (Fig.  3 A,  B),  the  least- 
squares  fits  fell  closer  to  the  1:1  line  than  the  eye  might  predict  beeausc  there  were  so  many 
overlapping  points.  At  412  nm  there  were  22,185  matched  points  and  at  526  nm  there  were 
28,491 ;  other  bands  fell  within  this  range.  Aggregate  plots  such  as  Fig.  3  gave  an  overview  of 
the  results,  but  details  of  the  models’  performance  were  apparent  when  looking  at  subsets  of 
the  data,  which  we  did  using  box  plots. 

Box  plots  show  the  distribution  of  model-data  differences,  D{0„  $,  across  the  observed 
chlorophyll  levels  (Figs.  4,  5),  azimuths  (Figs.  6,  7),  and  solar  zenith  angles  (Figs.  8,  9).  The 
plots  were  prepared  for  all  bands  but  are  only  shown  here  for  the  526  nm  channel  because  the 
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general  trends  were  the  same  at  all  wavelengths  (see  Fig.  3).  Box  plots  depict  the  distribution 
of  model-data  differences  as  follows:  the  box  shows  the  interquartile  range  (the  middle  50% 
of  the  data);  the  circle  with  dot  shows  the  median  model-data  difference;  thin  lines  known  as 
“whiskers”  extending  out  of  the  box  show  the  full  range  of  data  values.  Some  of  the 
distributions  of  D(0^  $  were  highly  skewed  or  contained  extreme  outliers;  box  plots  were 
useful  for  visualizing  such  features. 

The  population  (N)  of  each  distribution  is  given  along  the  top  of  each  plot.  Due  to  the  data 
acquired  with  chlorophyll  values  greater  than  10  pg  r\  the  number  of  model-data  differences 
was  not  always  the  same  for  the  different  models.  When  two  values  of  N  are  reported,  the  top 
one  refers  to  the  number  of  model-data  differences  for  the  MAG2002  look-up  table  and  the 
bottom  one  refers  to  the  number  of  model-data  differences  for  the  Lec20l  I  tables  and  our 
three  RT  model  calculations.  If  only  one  value  of  N  is  reported,  the  number  of  model-data 
differences  was  the  same  for  all  five  models. 

An  area  around  D(0„  <j>)  *  0  in  the  box  plots  has  been  shaded  to  show  +/-  the  mean 
coefficient  of  variation  of  the  NuRADS  images  at  the  points  used  to  compute  D(0„  <f>). 
Section  2.1  describes  how  the  CV  of  the  NuRADS  images  was  calculated.  The  CV  gives  a 
measure  of  the  limit  to  which  the  data  could  agree  with  even  the  best  model,  given 
environmental  variability,  due  to  effects  such  as  wave  focusing  within  the  images. 
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Rg.  3.  U0„  #)  /  U0,  0)  for  MAG2002  vs.  NuRADS  data  (A  and  C)  and  Lec20l 1  vs. 

NuRADS  data  (B  and  D).  The  plots  for  each  of  the  5  wavelengths  have  been  offset  vertically 
by  1 .0  for  clarity.  Solid  lines  on  each  plot  show  the  1 .1  slope  and  dashed  lines  show  a  linear  fit 
to  the  data.  Note  that  there  is  no  obvious  spectral  dependence  of  these  results. 

For  chlorophyll  concentrations  <  1 .0  pg  I  \  MAG2002  had  both  the  smallest  overall  range 
of  D{0„  i.e.  the  least  extreme  values,  and  the  smallest  interquartile  range  (Fig.  4). 
Furthermore,  over  this  range  MAG2002  generally  had  median  values  closest  to  zero,  although 
the  median  D{0„  $  for  Lcc20l  I  was  as  close  to  zero  as  MAG2002  for  chlorophyll 
concentrations  <  0. 1 8  pg  I"1  (Fig.  4). 


*  1 60488  -  SI  5.00  USD  Received  23  Dec  201 1 ;  revised  5  Mar  201 2;  accepted  5  Mar  2012;  published  20  Mar  201 2 
(C)20I2  OSA  26  March  2012/  Vol.  20,  No.  7  /  OPTICS  EXPRESS  7639 


At  chlorophyll  concentrations  >  1 .0  jag  I  1  in  the  Case  1  waters  we  sampled,  all  of  the 
models  had  similar  distributions  of  D{0y,  p),  It  was  not  obvious  that  one  model  performed 
consistently  better  than  the  others  for  chlorophyll  concentrations  >  1 .0  pg  I  1  (Fig.  4). 

Likewise,  in  the  samples  taken  from  Case  II  waters,  it  was  not  obvious  that  one  model 
performed  consistently  better  than  the  others  (Fig.  5).  In  the  chlorophyll  bins  for  which 
relatively  large  sample  sizes  were  availahle,  namely  1.78  -  3.15,  5.62  *  10.0,  10.0  -  17.8,  and 
31.6  -  56.2  jig  I  \  all  of  the  models  had  close  to  zero  median  and  similar  extreme  values  of 
D(0„  <f>).  In  the  chlorophyll  bins  with  smaller  sample  sizes,  1 .00  *  1 .78,  1 7.8  -  3 1 .6,  56.2  -  1 78, 
and  >  178  pg  I  \  all  of  the  models  exhibited  similar  biases  and  extreme  values  of 
approximately  equal  magnitude.  The  only  exceptions  to  these  generalizations  were  MAG2002 
and,  arguahly,  Lee201 1,  in  the  5.62  *  10.0  pg  I  1  bin.  These  differences,  however,  were  not 
part  of  any  trend.  MAG2002  and  Lec20l  I  performed  just  as  well  as  the  others  in  the  bin 
below  (1.78  -  3.16  pg  I  ')  and  the  bin  above  (10.0  -  17.8  pg  I  *)  the  one  in  question. 

In  Case  1  conditions,  all  of  the  models  had  positive  biases  in  and  largest  extreme  values  of 
D{6 U  p)  at  small  azimuth  angles  and  minimum  model-data  differences  between  105°  120° 

azimuth  (Fig.  6).  At  larger  azimuth  angles,  from  135°  to  180°,  the  model-data  agreement  for 
MAG2002  and  Lee201 1  did  not  decrease  much  relative  to  their  optimum  angles  between  105° 
-  120°  degrees,  but  the  biases  and  largest  model-data  deviations  for  the  other  models  steadily 
increased  (Fig.  6). 

MAG2002  had  the  smallest  bias  and  smallest  extreme  values  over  the  widest  range  of 
azimuth  angles  of  any  of  the  models  tested  (black  boxes  in  Fig.  6).  Our  RT  model  using  either 
the  Sullivan-Twardowski  (red  boxes)  or  Petzold  (blue  boxes)  phase  functions  scaled  to  give 
correct  bhp  values  had  approximately  the  same  bias,  though  slightly  larger  extreme  values,  as 
MAG2002  for  small  azimuth  angles,  but  they  had  notably  larger  biases  than  MAG2002  at 
large  azimuth  angles  (Fig.  6).  Conversely,  Lce20l  I  had  a  larger  bias  than  MAG2002  at  small 
azimuth  angles,  but  the  two  models  had  very  similar  performance  at  azimuth  angles  >  =  90°. 
Our  RT  model  using  the  Petzold  phase  function  scaled  to  give  correct  bp  values  had  larger 
biases  and  extreme  values  than  MAG2002  for  both  small  and  large  azimuth  angles  (green 
boxes  in  Fig.  6). 

One  way  in  which  the  Case  II  data  differed  from  the  Case  1  results,  however,  was  that 
MAG2002  exhibited  the  largest  biases  of  all  the  models,  whieh  were  found  from  0°  -45° 
azimuth.  A  second  way  in  which  the  Case  II  data  differed  from  the  Case  I  results  was  that  the 
range  of  azimuth  angles  over  which  there  was  good  model-data  agreement  was  much  larger 
for  the  Case  II  data.  For  example,  hoth  Lee20ll  and  our  RT  model  using  the  Sullivan- 
Twardowski  phase  function  produced  median  values  of  D(0„  p)  that  were  within  the  average 
NuRADS  coefficient  of  variation  across  all  azimuth  angles  (Fig.  7).  Third,  for  the  Case  II 
data,  the  range  of  extreme  values  of  D(0„  p)  fell  between  -10%  to  +  20%  of  the  upwelling 
radiance  at  nadir  for  azimuth  angles  45°  to  180°,  whereas  in  Case  I  the  smallest  extreme 
values  were  about  +/-  20%  only  between  azimuth  angles  of  90°  to  135°.  Finally,  between  0 
and  45°  the  largest  differences  were  positive  in  Case  I  and  negative  in  Case  II. 

Large  negative  biases  were  apparent  in  all  models  at  small  solar  zenith  angles  because  the 
NuRADS  data  were  normalized  to  nadir  (Fig.  8).  At  small  solar  zenith  angles,  near-nadir 
measurements  are  reduced  by  the  instrument  shadow  thereby  artificially  increasing  values  at 
other  angles  when  the  data  were  normalized  to  nadir.  For  the  Case  I  data,  the  near-nadir 
shadow  bias  was  eliminated  for  MAG2002  and  Lce201 1  for  0a  >  =  20°  and  for  our  RT  model 
for  0S  >  =  25°.  All  models  most  closely  agreed  with  the  data  at  solar  zenith  angles  between 
25°  to  35°  (Fig.  8).  At  0a  >  =  40°  a  positive  bias  in  D(0„  p)  and  highly  skewed  distributions 
were  observed  for  all  models  (Fig.  8).  The  median  value  of  D{0„  p)  computed  with 
MAG2002  exhibited  the  least  bias  for  0S  >  -  40°,  falling  within  the  average  NuRADS 
coefficient  of  variation  for  all  solar  zenith  angles  between  15°  to  70°.  Lce201 1  and  our  RT 
model  using  the  Sullivan-Twardowski  phase  function  produced  median  values  of  D{0„  P)  that 
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were  within  the  average  NuRADS  coefficient  of  variation  for  solar  zenith  angles  as  large  as 
65°. 


Fig  4.  Summary  of  1X0, ,  #)  as  a  function  of  chlorophyll  (Chi)  for  the  526  nm  NuRADS 
images  in  Case  1  water.  Sec  text  for  a  description  of  box  plots  and  of  the  population  sizes  (jV). 
From  left  to  right  within  each  group  of  boxes:  Black  boxes  correspond  to  values  from 
MAG2002  minus  the  data.  Light  blue  boxes  correspond  to  values  from  Lcc201 1  minus  the 
data.  Green,  red,  and  dark  blue  boxes  correspond  to  the  RT  model  using  VSF  parameterization 
methods  1-3,  respectively,  minus  the  data.  The  shaded  grey  area  around  0  difference  U  the 
mean  coefficient  of  variation  of  the  NuRADS  images  at  the  points  used  to  compute  1X0*  $ 


Fig.  5.  Summary  of  1X0*  #)  as  a  function  of  chlorophyll  (Chi)  for  the  526  nm  NuRADS 
images  in  Case  II  water.  See  text  and  caption  for  Fig.  4  for  description  of  boxes. 

The  Case  II  results,  plotted  against  solar  zenith  angle  (Fig.  9),  show  much  greater  model- 
data  agreement  as  well  as  consistency  among  the  models  than  the  data  from  Case  I  conditions 
(Ftg,  8).  The  only  examples  of  median  D(0V ,  </>)  falling  outside  the  range  of  the  average 


#160488  -  $1 5.00  USD  Received  23  Dec  201 1;  revised  5  Mar  2012;  accepted  5  Mar  2012;  published  20  Mar  2012 

(C)  20 1 2  OS  A  26  March  2012/  Vol.  20,  No.  7  /  OPTICS  EXPRESS  764 1 


NuRADS  coefficient  of  variation  were  for  situations  with  small  sample  sizes:  all  models  with 
35°  <  =  0,  <  40°  (N  =  29),  MAG2002  for  45°  <  =  0,  <  50 3  (N  =  9),  and  MAG2002  for  50°  <  = 
0S  <  55°  ( N  =  24).  No  systematic  patterns  of  bias  varying  with  solar  zenith  angle  were 
observed  (Fig.  9).  Note,  however,  that  there  were  no  observations  available  for  Case  II  water 
with  0t  <  35°.  Therefore  if  a  shadow  effect  such  as  the  one  observed  in  Case  I  water  for  0S  < 
25°  occurs  in  Case  II  water,  the  data  would  not  have  captured  it. 


526  nm  (Case  I) 

N  =  2322  2320  2322  2321  2320  2317  2317  2315  2300  2280  2194  1784  1379 


Fig.  6.  Summary  of  l)(0„  f)  as  a  function  of  view  azimuth  for  the  526  nm  NuRADS  images  in 
Case  1  water.  The  principal  plane  is  defined  by  azimuth  «  0°  (toward  the  sun)  and  azimuth  * 
180°  (away  from  the  sun).  See  text  and  caption  for  Rg.  4  for  description  of  boxes. 


526  nm  (Case  II) 


Rg.  7.  Summary  of  D(0„  ^)  as  a  function  of  view  azimuth  for  the  526  nm  NuRADS  images  in 
Case  II  water.  The  principal  plane  is  defined  by  azimuth  ■  0°  (toward  the  sun)  and  azimuth  * 
180°  (away  from  the  sun).  See  text  and  caption  for  Rg.  4  for  description  of  boxes. 
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Fig.  8.  Summary  of  l\0„  #)  as  a  function  of  0,  for  the  526  nm  NuRADS  images  in  Case  1 
water.  Sec  text  and  caption  for  Fig.  4  for  description  of  boxes. 
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Fig.  9.  Summary  of  D(0„  as  a  function  of  solar  zenith  angle  for  the  526  nm  NuRADS 
images  in  Case  11  water.  Sec  text  and  caption  for  Fig.  4  for  description  of  boxes. 

4.  Discussion 

Across  all  of  the  observed  Case  l  conditions,  MAG2002  agreed  with  much  of  the  available 
BRDF  data  more  closely  than  any  of  the  other  models.  In  particular,  three  portions  of  the  Case 
1  data  space  were  best  fit  by  MAG2002:  (a)  chlorophyll  concentrations  in  the  range  0.18  to 
1.0  pg  1  [  (Fig.  4),  (b)  azimuth  angles  near  the  principal  plane,  i.e.  outside  of  the  range  90^  - 
135°  (Fig.  6),  and  (e)  solar  zenith  angles  greater  than  35°  (Fig.  8).  At  both  low  and  high 
chlorophyll  concentrations,  and  for  solar  zenith  angles  from  -20-40  degrees,  however,  BRDF 
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correction  factors  from  Lcc201 1  and  our  method  2,  using  the  Sullivan  and  Twardowski  [13] 
phase  function,  each  matched  the  NuRADS  data  as  well  as  MAG2002.  All  of  the  models 
matched  the  NuRADS  data  well  for  azimuth  angles  near  90°.  It  is  not  surprising  that 
MAG2002  worked  well  in  Case  I  waters,  as  that  had  been  shown  previously  [II],  The  new 
results  here,  for  Case  I  water,  are  that  the  other,  models  with  simpler  parameterization  for 
their  particle  phase  functions  worked  nearly  as  well  as  MAG2002  under  some  conditions. 

For  Case  II  conditions,  in  contrast,  the  other  models  tested  matched  the  NuRADS  data  at 
least  as  well  as  the  corrections  from  MAG2002.  Lce20l  1  and  all  three  versions  of  our  RT 
code  were  hardly  distinguishable  in  Case  II  conditions  across  the  entire  range  of  chlorophyll 
concentrations  (Fig.  5),  azimuth  angles  (Fig.  7),  and  solar  zenith  angles  (Fig.  9)  sampled  by 
the  NuRADS  data.  The  RRDF  correction  generated  from  MAG2002  worked  well  in  Case  II 
water,  up  to  the  10  pg  I  1  limit  of  the  tables,  except  for  small  azimuth  angles  (Fig.  7).  Detailed 
results  have  been  presented  for  only  the  526  nm  data  (Figs.  4-9),  but  the  patterns  observed 
were  consistent  across  all  spectral  bands  sampled  (Fig.  3). 

Use  of  the  MAG2002  tables  without  Raman  scattering  would  not  likely  have  changed 
these  results  significantly.  MAG2002  results  with  and  without  Raman  scattering  differed  for / 
only  for  wavelengths  >  600  nm,  which  were  not  considered  in  this  study,  and  only  by  a  few 
percent  for  Q  [9]. 

The  agreement  of  our  RT  methods  1-3  with  the  data  revealed  that  both  the  total 
backscattering  (i.e.  the  value  of  b^)  and  the  detailed  shape  of  the  backward  portion  of  the 
particulate  phase  function  affected  the  BRDF  correction  in  Case  I  water.  At  small  azimuth 
angles,  methods  2  and  3,  in  which  the  value  of  bhp  matches  measurements,  both  agree  more 
closely  with  the  data  than  method  I,  in  which  the  value  of  bp  matches  measurements  (Fig.  6). 
This  result  indicates  that  the  magnitude  of  b^  is  more  important  for  BRDF  corrections  at 
small  azimuth  angles  than  the  detailed  shape  of  the  backward  portion  of  the  particulate  phase 
function. 

At  large  azimuth  angles,  however,  method  2,  using  the  Sullivan-Twardowski  phase 
function,  agreed  better  with  the  data  than  either  method  using  the  Petzold  phase  function  (Fig. 
6).  This  result  indicates  that  having  the  magnitude  of  b^  correct  was  necessary  but  not 
sufficient  to  match  the  data  at  large  azimuth  angles,  and  that  the  shape  of  the  particulate  phase 
function  in  the  backward  direction  was  also  important.  The  Sullivan-Twardowski  and  Petzold 
phase  functions  differ  by  about  5%  near  145°  and  109£  at  170°  (Fig.  I  inset). 

In  Case  II  waters,  the  differences  between  our  RT  methods  1-3  were  mueh  smaller  than  in 
Case  I  (Figs.  6,  7)  suggesting  that  multiple  scattering  dampened  out  the  most  of  the  effects  of 
the  phase  function  shape  for  determining  the  BRDF. 

The  observation  that  the  largest  model-data  discrepancies  were  close  to  the  principal  plane 
(Figs.  6,  7)  is  relevant  to  correcting  ocean  color  satellite  imagery  for  bidirectional  reflectance 
effects  because  the  area  of  largest  errors  is  also  the  portion  of  satellite  images  most  affected 
by  sunglint.  Ocean  color  processing  typically  discards  data  collected  within  and  close  to 
sunglint.  Therefore,  Moderate  Resolution  Imaging  Spcctroradiometer  (MODIS)  and  Sea- 
viewing  Wide  Field-of-view  Sensor  (SeaWiFS)  data  that  have  passed  quality  control  checks 
for  sunglint  will  have  been  acquired  over  azimuth  angles  with  the  smallest  model-data  BRDF 
discrepancies. 

5.  Conclusion 

Two  questions  addressed  were:  How  well  docs  the  recent  Lee 201 1  BRDF  correction  agree 
with  a  wide  suite  of  in  situ  measurements?  Is  the  choice  of  phase  function  used  in  Lcc201 1  a 
limitation  to  its  broad  implementation9  The  answers  to  these  questions  differed  in  Case  I  and 
Case  II  water. 

In  Case  I  conditions,  MAG2002  produced  the  closest  agreement  with  field  measurements 
of  the  models  considered  in  this  study.  Relative  to  earlier  work  [II]  that  also  validated  the 
accuracy  of  MAG2002  in  Case  I  water,  this  study  added  additional  data  sets  in  other  parts  of 
the  world  and  five  competing  models  for  comparison.  The  fact  that  MAG2002,  using  a 
variable  particulate  phase  function,  produced  the  best  results  and  that  the  other  models  with 
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different  particulate  phase  functions  did  not  always  produce  similar  output  suggests  that  in 
Case  1  water  the  shape  of  the  particulate  phase  function  is  important  for  modeling  BRDF  of 
the  ocean.  Specifically,  as  indicated  by  the  results  of  our  RT  methods  1-3,  the  detailed 
structure  of  the  backward  portion  of  the  particulate  phase  function  was  important  at  large 
azimuth  angles,  and  the  relative  amount  of  forward  to  backwards  scattering  was  important  at 
small  azimuth  angles. 

In  Case  II  conditions,  several  models  produced  similar  results;  of  the  ones  considered, 
Lce20l  1  is  probably  the  most  convenient  to  implement,  since  published  tables  are  available 
already,  Lce201 1  compared  their  model  with  two  NuRADS  images,  but  this  study  validated 
the  model  with  more  than  1000  NuRADS  images  across  a  wide  range  of  absorption  and 
scattering  coefficients,  solar  zenith  angles,  and  wavelengths.  Models  with  different,  but  fixed, 
phase  functions  produced  similar  output,  all  in  agreement  with  the  data,  suggesting  that  in 
Case  ii  water  the  precise  shape  of  the  phase  function  is  not  important  for  modeling  BRDF  of 
the  ocean  and  therefore  is  also  not  a  limitation  to  the  adoption  of  Lcc201 1  in  Case  li 
conditions. 

The  results  of  this  study  suggest  that  operational  correction  of  remotely  sensed  data  for 
bidirectional  reflectance  effects  could  be  performed  with  two  models,  depending  on  the  water 
type.  Lee201 1  (10)  was  superior  to  MAG2002  [9]  in  Case  II  water  at  small  azimuth  angles 
and  for  chlorophyll  concentrations  greater  than  the  10  pg  I  1  limit  of  MAG2002.  MAG2002 
[9],  however,  remained  superior  to  Lce20l  I  [10)  for  Case  I  water. 
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